####Figure Kd modeled vs Kd observed
setwd('/media/Iomega_HDD/boulot/Article/kd/last_23_01_2012/data/')

kd_lim<-6
kw<-c(0.030150000, 0.018650000, 0.009450000, 0.007632472, 0.009321161, 0.016500001, 0.033734146, 0.060471872)
wave=c(320,340,380,412,443,490,510,555)
nwave<-length(wave)

load("dataset.Rdata")
table <- cbind(signif(t5$value,7), signif(t5$gsm_chl,7), signif(t5$gsm_cdm,7),
                   signif(t5$gsm_bbp,7),t5$lambda, t5$convert_lambda)

para<-read.table('../script_admb/model.std',header=T, fill=T)
para <- para[,1:4]

s0<-para[9,3]
aphstar<-para[27:34,3]
a2<-para[18:25,3]
eta<-para[26,3]


dl<-matrix(0,nrow=1,ncol=3)
for (i in 1:nwave){
  kds<-table[which(table[,6]==i),1]		
	chl_gsm<-table[which(table[,6]==i),2]	
	acdm_gsm<-table[which(table[,6]==i),3]
	bbp_gsm<-table[which(table[,6]==i),4]	
	S<-s0
	CC<-acdm_gsm
	BBP<-bbp_gsm
	
	Kdmod<-kw[i]+CC*exp(-S*(wave[i]-443))+aphstar[i]*chl_gsm^a2[i]+BBP*(wave[i]/443)^(eta)
	titi<-cbind(kds,Kdmod, wave[i])
	dl<-rbind(dl,titi)
}
dl<-dl[2:length(dl[,1]),]


####Mean error

meanerr<-rep(0,8)
for (i in 1:8){
	
	toto<-subset(dl, subset=c(dl[,3]==wave[i]))
	tata<-abs((toto[,2]-toto[,1]))/toto[,1]*100
	tutu<-mean(tata)
	meanerr[i]<-tutu
	
	}


texti<-c("(a)", "(b)", "(c)", "(d)", "(e)", "(f)", "(g)", "(h)")
postscript('../Figures/Figure3.ps', width=10,height=15,paper="special", horizontal=F)
titick<-c(seq(0.01,0.1, by=0.01),seq(0.1,1, by=0.1), seq(1,5,by=1))
opar<-par()
par(mar=c(6,5,2,0.5))
par(mfrow=c(nwave/2,2))

for (i in 1:nwave){
plot(dl[which(dl[,3]==wave[i]),2]~dl[which(dl[,3]==wave[i]),1], pch=19,ylim=c(0.01,kd_lim),xlim=c(0.01,kd_lim), cex=0.7, cex.axis=2, cex.lab=2, xlab=expression(paste("Kd in situ", "  (",m^-1,")")), ylab=expression(paste("Kd predicted", "  (",m^-1,")")), main=paste(wave[i], " nm",sep=""),log="xy", xaxt="n", yaxt="n", cex.main=2)#
text(5,0.02, labels=texti[i], cex=2)
axis(side=1, at=titick, tick=T, labels=F)
axis(side=1, at=c(0.01,0.1,1,5), labels=T, tick=F, cex.axis=2)
axis(side=2, at=c(0.01,0.1,1,5), labels=T, tick=F, cex.axis=2)
axis(side=2, at=titick, tick=T, labels=F)
model<-lm(log(dl[which(dl[,3]==wave[i]),2])~log(dl[which(dl[,3]==wave[i]),1]))
model2<-lm(dl[which(dl[,3]==wave[i]),2]~dl[which(dl[,3]==wave[i]),1])
pupu<-summary(model)
print(summary(model2))
abline(a=0,b=1,lwd=2)
curve(x*0.5,0.001,kd_lim, col="dark green", lty=2, lwd=2, add=T)
curve(x*1.5,0.001,kd_lim, col="dark green", lty=2, lwd=2, add=T)
text(0.01,(kd_lim-0.1*kd_lim),labels=(paste("R²=",round(pupu$r.squared,digits=3),sep="")),cex=2,pos=4)
text(0.01,(kd_lim-0.5*kd_lim), labels=paste("slope=",round(model2$coefficients[2],digits=3),sep=""),cex=2,pos=4)
text(0.01,(kd_lim-0.75*kd_lim), labels=paste("N=",length(dl[which(dl[,3]==wave[i]),2]),sep=""),cex=2,pos=4) 
text(0.01,(kd_lim-0.87*kd_lim), labels=paste("mean error=",round(meanerr[i], digits=1),"%",sep=""),cex=2,pos=4)


}
dev.off()
par(opar)
